E16_Dec_subset

#load("E16_Dec7v3_Trachea.RData")

basal, secretory, and ciliated:

table(E16_Dec7v3_epi@meta.data$res.1.2)

   1   13   17   19   20    4    6    7    9 
1849  566  333  263  235 1362 1072  887  749 

colnames(E16_Dec7v3_epi@meta.data)[colnames(E16_Dec7v3_epi@meta.data) == 'res.0.8'] <- 'orig.0.8'
colnames(E16_Dec7v3_epi@meta.data)[colnames(E16_Dec7v3_epi@meta.data) == 'res.1.2'] <- 'orig.1.2'

E16_Dec7v3_epi <- ScaleData(object = E16_Dec7v3_epi)

Scaling data matrix

  |                                                                                                                                                               
  |                                                                                                                                                         |   0%
  |                                                                                                                                                               
  |=======                                                                                                                                                  |   5%
  |                                                                                                                                                               
  |==============                                                                                                                                           |   9%
  |                                                                                                                                                               
  |=====================                                                                                                                                    |  14%
  |                                                                                                                                                               
  |============================                                                                                                                             |  18%
  |                                                                                                                                                               
  |===================================                                                                                                                      |  23%
  |                                                                                                                                                               
  |==========================================                                                                                                               |  27%
  |                                                                                                                                                               
  |=================================================                                                                                                        |  32%
  |                                                                                                                                                               
  |========================================================                                                                                                 |  36%
  |                                                                                                                                                               
  |===============================================================                                                                                          |  41%
  |                                                                                                                                                               
  |======================================================================                                                                                   |  45%
  |                                                                                                                                                               
  |============================================================================                                                                             |  50%
  |                                                                                                                                                               
  |===================================================================================                                                                      |  55%
  |                                                                                                                                                               
  |==========================================================================================                                                               |  59%
  |                                                                                                                                                               
  |=================================================================================================                                                        |  64%
  |                                                                                                                                                               
  |========================================================================================================                                                 |  68%
  |                                                                                                                                                               
  |===============================================================================================================                                          |  73%
  |                                                                                                                                                               
  |======================================================================================================================                                   |  77%
  |                                                                                                                                                               
  |=============================================================================================================================                            |  82%
  |                                                                                                                                                               
  |====================================================================================================================================                     |  86%
  |                                                                                                                                                               
  |===========================================================================================================================================              |  91%
  |                                                                                                                                                               
  |==================================================================================================================================================       |  95%
  |                                                                                                                                                               
  |=========================================================================================================================================================| 100%

E16_Dec7v3_epi <- FindVariableGenes(object = E16_Dec7v3_epi, do.plot = TRUE, x.low.cutoff=0.1,x.high.cutoff = Inf, y.cutoff = 0.5)

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

run PCA on the set of genes

E16_Dec7v3_epi <- RunPCA(object = E16_Dec7v3_epi, do.print = FALSE)
#PCAPlot(E16_Dec7v3_epi)

E16_Dec7v3_epi <- ProjectPCA(object = E16_Dec7v3_epi, do.print = F)

n.pcs.sub = 17

res.used <- 1.2

E16_Dec7v3_epi <- FindClusters(object = E16_Dec7v3_epi, reduction.type = "pca", dims.use = 1:n.pcs.sub, 
                     resolution = res.used, print.output = 0, force.recalc = T)

E16_Dec7v3_epi <- RunTSNE(object = E16_Dec7v3_epi, dims.use = 1:n.pcs.sub, perplexity=30)

get tSNE embedding for velocyto:

E16_Dec_cv3_ID<-read.csv(file="E16_Dec_cv3_ID.csv",header=F,sep=",",stringsAsFactors = F) 

head(E16_Dec_cv3_ID)

E16_Dec_cv3_name<-gsub("x","",E16_Dec_cv3_ID)
head(E16_Dec_cv3_name)

E16_Dec_cv3_name<-gsub(":","_",E16_Dec_cv3_name)
head(E16_Dec_cv3_name)

TSNE1_Loomorder_epi<-E16_Dec7v3_epi@dr$tsne@cell.embeddings[match(E16_Dec_cv3_name,rownames(E16_Dec7v3_epi@dr$tsne@cell.embeddings)),1]
write(TSNE1_Loomorder_epi,"TSNE1_Loomorder_epi.csv",ncolumns=1,sep=",")

head(TSNE1_Loomorder_epi)

TSNE2_Loomorder_epi<-E16_Dec7v3_epi@dr$tsne@cell.embeddings[match(E16_Dec_cv3_name,rownames(E16_Dec7v3_epi@dr$tsne@cell.embeddings)),2]
write(TSNE2_Loomorder_epi,"TSNE2_Loomorder_epi.csv",ncolumns=1,sep=",")

markers for each ciliated population:

E16_Dec_epi_res12_8over1011<-FindMarkers(E16_Dec7v3_epi,ident.1=c(8),ident.2 = c(10,11),only.pos = TRUE)
E16_Dec_epi_res12_8over1011

write.table(E16_Dec_epi_res12_8over1011,"epiSubset_c8inCilia.txt",sep="\t")

E16_Dec_epi_res12_10over811<-FindMarkers(E16_Dec7v3_epi,ident.1=c(10),ident.2 = c(8,11),only.pos = TRUE)
E16_Dec_epi_res12_10over811

E16_Dec_epi_res12_10over8<-FindMarkers(E16_Dec7v3_epi,ident.1=c(10),ident.2 = c(8),only.pos = TRUE)
E16_Dec_epi_res12_10over8

write.table(E16_Dec_epi_res12_10over8,"epiSubset_c10overC8.txt",sep="\t")

E16_Dec7v3_epi <- SetAllIdent(object = E16_Dec7v3_epi, id = "res.1.2")

E16_Dec_epi_res12_10over11<-FindMarkers(E16_Dec7v3_epi,ident.1=c(10),ident.2 = c(11),only.pos = TRUE)
E16_Dec_epi_res12_10over11

write.table(E16_Dec_epi_res12_10over11,"epiSubset_c10overC11.txt",sep="\t")

E16_Dec_epi_res12_11over810<-FindMarkers(E16_Dec7v3_epi,ident.1=c(11),ident.2 = c(8,10),only.pos = TRUE)
E16_Dec_epi_res12_11over810

write.table(E16_Dec_epi_res12_11over810,"epiSubset_c11inCilia.txt",sep="\t")

E16_Dec7v3_epi <- SetAllIdent(object = E16_Dec7v3_epi, id = "res.1.2")
DoHeatmap(object = E16_Dec7v3_epi, genes.use = c("Foxj1","Top2a","Mcidas","Ccno","Foxn4","Shisa8","Lrrc23","Prr18","Cfap53","Cdhr3","Sntn","Ifitm1","Lbp","Ly6c1","Ly6a"), 
    slim.col.label = TRUE, group.label.rot = TRUE,use.scaled = F,cells.use = E16_Dec7v3_epi@cell.names[E16_Dec7v3_epi@meta.data$res.1.2 %in% c(8,10,11)],group.order = c(8,10,11),group.cex = 30,cex.row = 20
  )

scoring:

E16_Dec7v3_epi@data[1:6,1:6]

percentile_table_epi<-apply(E16_Dec7v3_epi@data,1,percent_rank)

percentile_table_epi[1:6,1:6]

OMIMgene<-read.csv(file = "genesOMIM.csv",header=T,sep=",",stringsAsFactors = F)
OMIMgene<-lapply(OMIMgene,function(x) unlist(strsplit(unlist(x),split=","))) 
head(OMIMgene$Mucociliary)

OMIMgene_mucosaGoblet<-as.vector(read.csv(file = "genesOMIM_mucosa_goblet.csv",header=T,sep=",",stringsAsFactors = F)[,1])
OMIMgene_mucosaGoblet<-unlist(strsplit(unlist(OMIMgene_mucosaGoblet),split=","))
OMIMgene_mucosaGoblet[90:105]

 mocosaGoblet_score<- apply(percentile_table_epi[,colnames(percentile_table_epi) %in% OMIMgene_mucosaGoblet],1,mean)

 head( mocosaGoblet_score)

E16_Dec7v3_epi<-AddMetaData(object = E16_Dec7v3_epi, metadata = mocosaGoblet_score, col.name = "mocosaGoblet_score")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("mocosaGoblet_score"), nCol = 1,x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="res.1.4")

 ciliopathy_table<- percentile_table_epi[,colnames(percentile_table_epi) %in% OMIMgene$Ciliopathy]

 ciliopathy_score<- apply(ciliopathy_table,1,mean)

 head(ciliopathy_score)

 PCD_score<- apply(percentile_table_epi[,colnames(percentile_table_epi) %in% OMIMgene$Primary.ciliary.dyskinesia],1,mean)

 head(PCD_score)

E16_Dec7v3_epi<-AddMetaData(object = E16_Dec7v3_epi, metadata = ciliopathy_score, col.name = "ciliopathy_score")

E16_Dec7v3_epi<-AddMetaData(object = E16_Dec7v3_epi, metadata = PCD_score, col.name = "PCD_score")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("ciliopathy_score"), nCol = 1,x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="res.1.2")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("PCD_score"), nCol = 1,x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="res.1.2")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("PCD_score"), nCol = 1,ident.include = c(8,10,11),x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="res.1.2")

E16_Dec7v3_epi <- SetAllIdent(object = E16_Dec7v3_epi, id = "res.1.2")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("ciliopathy_score"), ident.include = c(11),nCol = 1,x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="seq_group")

E16_Dec7v3_epi <- SetAllIdent(object = E16_Dec7v3_epi, id = "res.1.2")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("PCD_score"), ident.include = c(8),nCol = 1,x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="seq_group")

 mucus_score<- apply(percentile_table_epi[,colnames(percentile_table_epi) %in% c(OMIMgene$Airway...Mucus,OMIMgene$Pulmonary.and.Mucus)],1,mean)

 head(mucus_score)

E16_Dec7v3_epi<-AddMetaData(object = E16_Dec7v3_epi, metadata = mucus_score, col.name = "mucus_score")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("mucus_score"), nCol = 1,x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="res.1.4")

 COPD_score<- apply(percentile_table_epi[,colnames(percentile_table_epi) %in% c(OMIMgene$COPD)],1,mean)

 head(COPD_score)

E16_Dec7v3_epi<-AddMetaData(object = E16_Dec7v3_epi, metadata = COPD_score, col.name = "COPD_score")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("COPD_score"), nCol = 1,x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="res.1.4")

 asthma_score<- apply(percentile_table_epi[,colnames(percentile_table_epi) %in% c(OMIMgene$Pulmonary...Asthma)],1,mean)

 head(asthma_score)

E16_Dec7v3_epi<-AddMetaData(object = E16_Dec7v3_epi, metadata = asthma_score, col.name = "asthma_score")

VlnPlot(object = E16_Dec7v3_epi, features.plot = c("asthma_score"), nCol = 1,x.lab.rot = T,point.size.use = 0.3,use.raw=F,group.by="res.1.4")

ggplot(E16_Dec7v3_epi@meta.data,aes(genotype,asthma_score))+facet_grid(.~res.1.2)+geom_dotplot(binaxis="y",aes(color=genotype,fill=genotype),binwidth=0.05,stackdir="center",position=position_dodge(0.8), dotsize=0.2)+stat_compare_means(comparisons = list(c("wt", "mut")),method="wilcox.test",size=4,label="p.adj")+ stat_summary(aes(color=genotype),fun.data=mean_sdl, fun.args = list(mult=1), 
                 geom="pointrange",position=position_dodge(0.7))+ theme(axis.text.x = element_text(angle = 45,hjust=1))

res.used <- 1.4

E16_Dec7v3_epi <- FindClusters(object = E16_Dec7v3_epi, reduction.type = "pca", dims.use = 1:n.pcs.sub, 
                     resolution = res.used, print.output = 0, force.recalc = T)

E16_Dec7v3_epi <- RunTSNE(object = E16_Dec7v3_epi, dims.use = 1:n.pcs.sub, perplexity=30)

E16_Dec7v3_epi=buildClusterTree(E16_Dec7v3_epi,do.reorder = F,reorder.numeric = F,pcs.use = 1:17)

table(E16_Dec7v3_epi@meta.data$res.1.4,E16_Dec7v3_epi@meta.data$seq_group)

prop.table(table(E16_Dec7v3_epi@meta.data$res.1.4,E16_Dec7v3_epi@meta.data$seq_group),2)

    
     E16_Dec7_mut_7 E16_Dec7_mut_8 E16_Dec7_wt_1 E16_Dec7_wt_6
  0     0.032475749    0.059523810   0.242384964   0.153679654
  1     0.113454239    0.068948413   0.099157485   0.145743146
  10    0.021931674    0.095238095   0.032404407   0.014430014
  11    0.050611556    0.028273810   0.036292936   0.048340548
  12    0.043441586    0.033730159   0.036292936   0.041125541
  13    0.032053986    0.026289683   0.043421905   0.028860029
  14    0.024040489    0.039682540   0.014906027   0.038239538
  15    0.051455082    0.022321429   0.005832793   0.014430014
  16    0.009700548    0.041666667   0.011665587   0.041847042
  17    0.029523408    0.024305556   0.002592353   0.020923521
  18    0.012652889    0.011904762   0.009073234   0.012987013
  2     0.099114298    0.112599206   0.050550875   0.084415584
  3     0.074230283    0.051091270   0.106934543   0.093073593
  4     0.088570224    0.096726190   0.053143227   0.057720058
  5     0.043019823    0.086805556   0.079066753   0.059884560
  6     0.147195276    0.023809524   0.030460143   0.023809524
  7     0.060733868    0.041666667   0.059624109   0.043290043
  8     0.014339941    0.082837302   0.053791316   0.046176046
  9     0.051455082    0.052579365   0.032404407   0.031024531

E16_Dec7v3_epi@meta.data$cell_type<-mapvalues(E16_Dec7v3_epi@meta.data$res.1.4,from=c("0","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18"),to=c("Basal","Basal","Secretory","Basal/Secretory","Secretory","Secretory","Secretory","Ciliated","Secretory","Secretory","Secretory","Ciliated","Secretory","Basal","Doublet","Ciliated","Secretory","Doublet","Doublet"))

c14, 17 and c18 are doublets.

table(E16_Dec7v3_epi@meta.data$cell_type[!(E16_Dec7v3_epi@meta.data$res.1.4 %in% c(14,17,18))],E16_Dec7v3_epi@meta.data$genotype[!(E16_Dec7v3_epi@meta.data$res.1.4 %in% c(14,17,18))])

                 
                   mut   wt
  Basal            734 1049
  Basal/Secretory  279  294
  Ciliated         572  304
  Secretory       2492 1141

DE_E16_secretory_genotype<-FindMarkers(E16_Dec7v3_epi,cells.1<-WhichCells(object=E16_Dec7v3_epi,cells.use = (E16_Dec7v3_epi@meta.data$genotype=="wt" & E16_Dec7v3_epi@meta.data$cell_type=="Secretory" )),cells.2<-WhichCells(object=E16_Dec7v3_epi,cells.use = (E16_Dec7v3_epi@meta.data$genotype=="mut" & E16_Dec7v3_epi@meta.data$cell_type=="Secretory" )),only.pos = F,logfc.threshold=0,min.pct=0)

   |                                                  | 0 % ~calculating  
   |+                                                 | 1 % ~30m 46s      
   |+                                                 | 2 % ~30m 22s      
   |++                                                | 3 % ~30m 21s      
   |++                                                | 4 % ~30m 15s      
   |+++                                               | 5 % ~29m 42s      
   |+++                                               | 6 % ~29m 17s      
   |++++                                              | 7 % ~28m 51s      
   |++++                                              | 8 % ~28m 30s      
   |+++++                                             | 9 % ~28m 06s      
   |+++++                                             | 10% ~27m 47s      
   |++++++                                            | 11% ~27m 28s      
   |++++++                                            | 12% ~27m 08s      
   |+++++++                                           | 13% ~26m 51s      
   |+++++++                                           | 14% ~26m 32s      
   |++++++++                                          | 15% ~26m 14s      
   |++++++++                                          | 16% ~25m 55s      
   |+++++++++                                         | 17% ~25m 37s      
   |+++++++++                                         | 18% ~25m 17s      
   |++++++++++                                        | 19% ~24m 58s      
   |++++++++++                                        | 20% ~24m 43s      
   |+++++++++++                                       | 21% ~24m 25s      
   |+++++++++++                                       | 22% ~24m 07s      
   |++++++++++++                                      | 23% ~23m 51s      
   |++++++++++++                                      | 24% ~23m 31s      
   |+++++++++++++                                     | 25% ~23m 12s      
   |+++++++++++++                                     | 26% ~22m 52s      
   |++++++++++++++                                    | 27% ~22m 34s      
   |++++++++++++++                                    | 28% ~22m 15s      
   |+++++++++++++++                                   | 29% ~21m 55s      
   |+++++++++++++++                                   | 30% ~21m 37s      
   |++++++++++++++++                                  | 31% ~21m 18s      
   |++++++++++++++++                                  | 32% ~20m 59s      
   |+++++++++++++++++                                 | 33% ~20m 40s      
   |+++++++++++++++++                                | 34% ~20m 22s      
   |++++++++++++++++++                                | 35% ~20m 04s      
   |++++++++++++++++++                                | 36% ~19m 46s      
   |+++++++++++++++++++                               | 37% ~19m 28s      
   |+++++++++++++++++++                               | 38% ~19m 10s      
   |++++++++++++++++++++                              | 39% ~18m 51s      
   |++++++++++++++++++++                              | 40% ~18m 33s      
   |+++++++++++++++++++++                             | 41% ~18m 15s      
   |+++++++++++++++++++++                             | 42% ~17m 57s      
   |++++++++++++++++++++++                            | 43% ~17m 38s      
   |++++++++++++++++++++++                            | 44% ~17m 19s      
   |+++++++++++++++++++++++                           | 45% ~17m 01s      
   |+++++++++++++++++++++++                           | 46% ~16m 42s      
   |++++++++++++++++++++++++                          | 47% ~16m 24s      
   |++++++++++++++++++++++++                          | 48% ~16m 05s      
   |+++++++++++++++++++++++++                         | 49% ~15m 46s      
   |+++++++++++++++++++++++++                         | 50% ~15m 27s      
   |++++++++++++++++++++++++++                        | 51% ~15m 08s      
   |++++++++++++++++++++++++++                        | 52% ~14m 50s      
   |+++++++++++++++++++++++++++                       | 53% ~14m 31s      
   |+++++++++++++++++++++++++++                       | 54% ~14m 13s      
   |++++++++++++++++++++++++++++                      | 55% ~13m 54s      
   |++++++++++++++++++++++++++++                     | 56% ~13m 35s      
   |+++++++++++++++++++++++++++++                     | 57% ~13m 17s      
   |+++++++++++++++++++++++++++++                     | 58% ~12m 58s      
   |++++++++++++++++++++++++++++++                    | 59% ~12m 40s      
   |++++++++++++++++++++++++++++++                    | 60% ~12m 22s      
   |+++++++++++++++++++++++++++++++                   | 61% ~12m 03s      
   |+++++++++++++++++++++++++++++++                   | 62% ~11m 45s      
   |++++++++++++++++++++++++++++++++                  | 63% ~11m 26s      
   |++++++++++++++++++++++++++++++++                  | 64% ~11m 08s      
   |+++++++++++++++++++++++++++++++++                 | 65% ~10m 49s      
   |+++++++++++++++++++++++++++++++++                 | 66% ~10m 31s      
   |++++++++++++++++++++++++++++++++++                | 67% ~10m 12s      
   |++++++++++++++++++++++++++++++++++               | 68% ~09m 53s      
   |+++++++++++++++++++++++++++++++++++               | 69% ~09m 35s      
   |+++++++++++++++++++++++++++++++++++               | 70% ~09m 16s      
   |++++++++++++++++++++++++++++++++++++              | 71% ~08m 58s      
   |++++++++++++++++++++++++++++++++++++              | 72% ~08m 39s      
   |+++++++++++++++++++++++++++++++++++++             | 73% ~08m 20s      
   |+++++++++++++++++++++++++++++++++++++             | 74% ~08m 02s      
   |++++++++++++++++++++++++++++++++++++++            | 75% ~07m 43s      
   |++++++++++++++++++++++++++++++++++++++            | 76% ~07m 25s      
   |+++++++++++++++++++++++++++++++++++++++           | 77% ~07m 06s      
   |+++++++++++++++++++++++++++++++++++++++          | 78% ~06m 48s      
   |++++++++++++++++++++++++++++++++++++++++          | 79% ~06m 29s      
   |++++++++++++++++++++++++++++++++++++++++         | 80% ~06m 11s      
   |+++++++++++++++++++++++++++++++++++++++++         | 81% ~05m 53s      
   |+++++++++++++++++++++++++++++++++++++++++         | 82% ~05m 34s      
   |++++++++++++++++++++++++++++++++++++++++++        | 83% ~05m 15s      
   |++++++++++++++++++++++++++++++++++++++++++        | 84% ~04m 57s      
   |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~04m 38s      
   |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~04m 20s      
   |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~04m 01s      
   |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~03m 42s      
   |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~03m 24s      
   |+++++++++++++++++++++++++++++++++++++++++++++    | 90% ~03m 05s      
   |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~02m 47s      
   |++++++++++++++++++++++++++++++++++++++++++++++   | 92% ~02m 28s      
   |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~02m 10s      
   |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~01m 51s      
   |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~01m 33s      
   |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~01m 14s      
   |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~56s          
   |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~37s          
   |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~19s          
   |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 30m 50s

DE_E16_secretory_genotype

library(ggrepel)

DE_E16_secretory_genotype$gene<-rownames(DE_E16_secretory_genotype)
DE_E16_secretory_genotype$sig<-DE_E16_secretory_genotype$p_val_adj<0.001
volc = ggplot(DE_E16_secretory_genotype, aes(avg_logFC, -log10(p_val_adj))) + #volcanoplot with avg_logFC versus p_val_adj
    geom_point(aes(col=sig)) + #add points colored by significance
    scale_color_manual(values=c("black", "red")) + 
    ggtitle("E16secretory_wt/mut") + geom_text_repel(data=head(DE_E16_secretory_genotype, 20), aes(label=gene), point.padding = 1, box.padding = .3) +
  labs(y = expression(-log[10]*" "*"adjusted pvalue"), x = "avg log fold change") + 
  theme(legend.title = element_blank(), legend.position = "top") + 
  scale_fill_discrete(labels = c("Not Sig", "adjusted pval < 0.001"))

volc

E16_Dec7v3_epi@meta.data$specific_type<-mapvalues(E16_Dec7v3_epi@meta.data$res.1.4,from=c("0","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18"),to=c("Basal-Krt14","Basal","Secretory-Krt4","Basal/Secretory","Secretory-Krt4","Secretory","Secretory-Krt4","Ciliated","Secretory","Secretory-Krt4","Secretory-Krt4","Ciliated","CyclingSecretory","CyclingBasal","Doublet","Ciliated","Secretory","Doublet","Doublet"))

table(E16_Dec7v3_epi@meta.data$specific_type[!(E16_Dec7v3_epi@meta.data$res.1.4 %in% c(14,17,18))],E16_Dec7v3_epi@meta.data$genotype[!(E16_Dec7v3_epi@meta.data$res.1.4 %in% c(14,17,18))])

                  
                    mut   wt
  Basal             408  355
  Basal/Secretory   279  294
  Basal-Sostdc1     197  587
  Ciliated          572  304
  CyclingBasal      129  107
  CyclingSecretory  171  113
  Secretory         585  428
  Secretory-Krt4   1736  600

markers for specific clusters:

print(levels(E16_Dec7v3_epi@ident))

[1] "Basal"            "Basal/Secretory"  "Basal-Sostdc1"    "Ciliated"         "CyclingBasal"     "CyclingSecretory" "Doublet"          "Secretory"       
[9] "Secretory-Krt4"  

different populations of basal cells:

write.table(E16_Dec_epi_res14_0over1,"epiSubset_res14_c0overC1.txt",sep="\t")

library(ggrepel)

E16_Dec_epi_res14_0over1$gene<-rownames(E16_Dec_epi_res14_0over1)
E16_Dec_epi_res14_0over1$sig<-E16_Dec_epi_res14_0over1$p_val_adj<0.001
volc = ggplot(E16_Dec_epi_res14_0over1, aes(avg_logFC, -log10(p_val_adj))) + #volcanoplot with avg_logFC versus p_val_adj
    geom_point(aes(col=sig)) + #add points colored by significance
    scale_color_manual(values=c("black", "red")) + 
    ggtitle("E16_epi_0over1") 

E16_Dec_epi_res14_0over1_volc.1 <- volc + geom_text_repel(data=head(E16_Dec_epi_res14_0over1, 20), aes(label=gene), point.padding = 1, box.padding = .3) +
  labs(y = expression(-log[10]*" "*"adjusted pvalue"), x = "avg log fold change") + 
  theme(legend.title = element_blank(), legend.position = "top") + 
  scale_fill_discrete(labels = c("Not Sig", "adjusted pval < 0.001"))

E16_Dec7v3_epi <- SetAllIdent(object = E16_Dec7v3_epi, id = "res.1.4")
E16_Dec_epi_res14_1over0<-FindMarkers(E16_Dec7v3_epi,ident.1=c(1),ident.2 = c(0),only.pos = TRUE)

   |                                                  | 0 % ~calculating  
   |+                                                 | 1 % ~05s          
   |++                                                | 3 % ~05s          
   |++                                                | 4 % ~04s          
   |+++                                               | 5 % ~03s          
   |++++                                              | 6 % ~03s          
   |++++                                              | 8 % ~03s          
   |+++++                                             | 9 % ~03s          
   |++++++                                            | 10% ~03s          
   |++++++                                            | 11% ~02s          
   |+++++++                                           | 13% ~02s          
   |+++++++                                           | 14% ~02s          
   |++++++++                                          | 15% ~02s          
   |+++++++++                                         | 16% ~02s          
   |+++++++++                                         | 18% ~02s          
   |++++++++++                                        | 19% ~02s          
   |+++++++++++                                       | 20% ~02s          
   |+++++++++++                                       | 22% ~02s          
   |++++++++++++                                      | 23% ~02s          
   |+++++++++++++                                     | 24% ~02s          
   |+++++++++++++                                     | 25% ~02s          
   |++++++++++++++                                    | 27% ~02s          
   |++++++++++++++                                    | 28% ~02s          
   |+++++++++++++++                                   | 29% ~02s          
   |++++++++++++++++                                  | 30% ~02s          
   |++++++++++++++++                                  | 32% ~02s          
   |+++++++++++++++++                                 | 33% ~02s          
   |++++++++++++++++++                                | 34% ~02s          
   |++++++++++++++++++                                | 35% ~02s          
   |+++++++++++++++++++                               | 37% ~01s          
   |+++++++++++++++++++                               | 38% ~01s          
   |++++++++++++++++++++                              | 39% ~01s          
   |+++++++++++++++++++++                             | 41% ~01s          
   |+++++++++++++++++++++                             | 42% ~01s          
   |++++++++++++++++++++++                            | 43% ~01s          
   |+++++++++++++++++++++++                           | 44% ~01s          
   |+++++++++++++++++++++++                           | 46% ~01s          
   |++++++++++++++++++++++++                          | 47% ~01s          
   |+++++++++++++++++++++++++                         | 48% ~01s          
   |+++++++++++++++++++++++++                         | 49% ~01s          
   |++++++++++++++++++++++++++                        | 51% ~01s          
   |++++++++++++++++++++++++++                        | 52% ~01s          
   |+++++++++++++++++++++++++++                       | 53% ~01s          
   |++++++++++++++++++++++++++++                      | 54% ~01s          
   |++++++++++++++++++++++++++++                      | 56% ~01s          
   |+++++++++++++++++++++++++++++                     | 57% ~01s          
   |++++++++++++++++++++++++++++++                    | 58% ~01s          
   |++++++++++++++++++++++++++++++                    | 59% ~01s          
   |+++++++++++++++++++++++++++++++                   | 61% ~01s          
   |++++++++++++++++++++++++++++++++                  | 62% ~01s          
   |++++++++++++++++++++++++++++++++                  | 63% ~01s          
   |+++++++++++++++++++++++++++++++++                 | 65% ~01s          
   |+++++++++++++++++++++++++++++++++                 | 66% ~01s          
   |++++++++++++++++++++++++++++++++++                | 67% ~01s          
   |+++++++++++++++++++++++++++++++++++               | 68% ~01s          
   |+++++++++++++++++++++++++++++++++++               | 70% ~01s          
   |++++++++++++++++++++++++++++++++++++              | 71% ~01s          
   |+++++++++++++++++++++++++++++++++++++             | 72% ~01s          
   |+++++++++++++++++++++++++++++++++++++             | 73% ~01s          
   |++++++++++++++++++++++++++++++++++++++            | 75% ~01s          
   |++++++++++++++++++++++++++++++++++++++            | 76% ~01s          
   |+++++++++++++++++++++++++++++++++++++++           | 77% ~01s          
   |++++++++++++++++++++++++++++++++++++++++          | 78% ~01s          
   |++++++++++++++++++++++++++++++++++++++++          | 80% ~00s          
   |+++++++++++++++++++++++++++++++++++++++++         | 81% ~00s          
   |++++++++++++++++++++++++++++++++++++++++++        | 82% ~00s          
   |++++++++++++++++++++++++++++++++++++++++++        | 84% ~00s          
   |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~00s          
   |++++++++++++++++++++++++++++++++++++++++++++      | 86% ~00s          
   |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~00s          
   |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~00s          
   |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~00s          
   |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~00s          
   |+++++++++++++++++++++++++++++++++++++++++++++++   | 92% ~00s          
   |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~00s          
   |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~00s          
   |+++++++++++++++++++++++++++++++++++++++++++++++++ | 96% ~00s          
   |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s          
   |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s          
   |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 02s

E16_Dec_epi_res14_1over0

write.table(E16_Dec_epi_res14_1over0,"epiSubset_res14_c1overC0.txt",sep="\t")

E16_Dec7v3_epi <- SetAllIdent(object = E16_Dec7v3_epi, id = "res.1.4")
E16_Dec_epi_res14_9over1<-FindMarkers(E16_Dec7v3_epi,ident.1=c(9),ident.2 = c(1),only.pos = TRUE)
E16_Dec_epi_res14_9over1

write.table(E16_Dec_epi_res14_9over1,"epiSubset_res14_c9overC1.txt",sep="\t")

E16_Dec7v3_epi <- SetAllIdent(object = E16_Dec7v3_epi, id = "res.1.4")
E16_Dec_epi_res14_1over9<-FindMarkers(E16_Dec7v3_epi,ident.1=c(1),ident.2 = c(9),only.pos = TRUE)
E16_Dec_epi_res14_1over9

write.table(E16_Dec_epi_res14_1over9,"epiSubset_res14_c1overC9.txt",sep="\t")

different populations of ciliated cells:

different populations of secretory cells:

Sostdc1 seems to be downregulated in mutant basal-Sostdc1 population.

epi_marker_c0_mut_wt<-FindMarkers(E16_Dec7v3_epi,cells.1<-WhichCells(object=E16_Dec7v3_epi,ident=0,cells.use = E16_Dec7v3_epi@meta.data$genotype=="mut"),cells.2<-WhichCells(object=E16_Dec7v3_epi,ident=0,cells.use = E16_Dec7v3_epi@meta.data$genotype=="wt"),only.pos = TRUE)
epi_marker_c0_mut_wt

E16_Dec7v3_epi@meta.data$type_genotype<-as.factor(paste(E16_Dec7v3_epi@meta.data$cell_type,E16_Dec7v3_epi@meta.data$genotype,sep="_"))

DE_E16_ciliated_genotype<-FindMarkers(E16_Dec7v3_epi,cells.1<-WhichCells(object=E16_Dec7v3_epi,cells.use = (E16_Dec7v3_epi@meta.data$genotype=="wt" & E16_Dec7v3_epi@meta.data$cell_type=="Ciliated" )),cells.2<-WhichCells(object=E16_Dec7v3_epi,cells.use = (E16_Dec7v3_epi@meta.data$genotype=="mut" & E16_Dec7v3_epi@meta.data$cell_type=="Ciliated" )),only.pos = F,logfc.threshold=0,min.pct=0.05)

   |                                                  | 0 % ~calculating  
   |+                                                 | 1 % ~05m 51s      
   |+                                                 | 2 % ~05m 39s      
   |++                                                | 3 % ~05m 36s      
   |++                                                | 4 % ~05m 29s      
   |+++                                               | 5 % ~05m 29s      
   |+++                                               | 6 % ~05m 23s      
   |++++                                              | 7 % ~05m 20s      
   |++++                                              | 8 % ~05m 15s      
   |+++++                                             | 9 % ~05m 14s      
   |+++++                                             | 10% ~05m 09s      
   |++++++                                            | 11% ~05m 05s      
   |++++++                                            | 12% ~05m 04s      
   |+++++++                                           | 13% ~05m 02s      
   |+++++++                                           | 14% ~04m 57s      
   |++++++++                                          | 15% ~04m 54s      
   |++++++++                                          | 16% ~04m 50s      
   |+++++++++                                         | 17% ~04m 47s      
   |+++++++++                                         | 18% ~04m 42s      
   |++++++++++                                        | 19% ~04m 39s      
   |++++++++++                                        | 20% ~04m 35s      
   |+++++++++++                                       | 21% ~04m 32s      
   |+++++++++++                                       | 22% ~04m 28s      
   |++++++++++++                                      | 23% ~04m 25s      
   |++++++++++++                                      | 24% ~04m 21s      
   |+++++++++++++                                     | 25% ~04m 18s      
   |+++++++++++++                                     | 26% ~04m 14s      
   |++++++++++++++                                    | 27% ~04m 11s      
   |++++++++++++++                                    | 28% ~04m 08s      
   |+++++++++++++++                                   | 29% ~04m 05s      
   |+++++++++++++++                                   | 30% ~04m 01s      
   |++++++++++++++++                                  | 31% ~03m 58s      
   |++++++++++++++++                                  | 32% ~03m 55s      
   |+++++++++++++++++                                 | 33% ~03m 51s      
   |+++++++++++++++++                                | 34% ~03m 47s      
   |++++++++++++++++++                                | 35% ~03m 44s      
   |++++++++++++++++++                                | 36% ~03m 40s      
   |+++++++++++++++++++                               | 37% ~03m 37s      
   |+++++++++++++++++++                               | 38% ~03m 34s      
   |++++++++++++++++++++                              | 39% ~03m 30s      
   |++++++++++++++++++++                              | 40% ~03m 27s      
   |+++++++++++++++++++++                             | 41% ~03m 24s      
   |+++++++++++++++++++++                             | 42% ~03m 20s      
   |++++++++++++++++++++++                            | 43% ~03m 18s      
   |++++++++++++++++++++++                            | 44% ~03m 14s      
   |+++++++++++++++++++++++                           | 45% ~03m 10s      
   |+++++++++++++++++++++++                           | 46% ~03m 07s      
   |++++++++++++++++++++++++                          | 47% ~03m 03s      
   |++++++++++++++++++++++++                          | 48% ~02m 59s      
   |+++++++++++++++++++++++++                         | 49% ~02m 56s      
   |+++++++++++++++++++++++++                         | 50% ~02m 52s      
   |++++++++++++++++++++++++++                        | 51% ~02m 49s      
   |++++++++++++++++++++++++++                        | 52% ~02m 45s      
   |+++++++++++++++++++++++++++                       | 53% ~02m 41s      
   |+++++++++++++++++++++++++++                       | 54% ~02m 38s      
   |++++++++++++++++++++++++++++                      | 55% ~02m 35s      
   |++++++++++++++++++++++++++++                     | 56% ~02m 31s      
   |+++++++++++++++++++++++++++++                     | 57% ~02m 28s      
   |+++++++++++++++++++++++++++++                     | 58% ~02m 24s      
   |++++++++++++++++++++++++++++++                    | 59% ~02m 21s      
   |++++++++++++++++++++++++++++++                    | 60% ~02m 17s      
   |+++++++++++++++++++++++++++++++                   | 61% ~02m 14s      
   |+++++++++++++++++++++++++++++++                   | 62% ~02m 10s      
   |++++++++++++++++++++++++++++++++                  | 63% ~02m 07s      
   |++++++++++++++++++++++++++++++++                  | 64% ~02m 03s      
   |+++++++++++++++++++++++++++++++++                 | 65% ~01m 60s      
   |+++++++++++++++++++++++++++++++++                 | 66% ~01m 56s      
   |++++++++++++++++++++++++++++++++++                | 67% ~01m 53s      
   |++++++++++++++++++++++++++++++++++               | 68% ~01m 49s      
   |+++++++++++++++++++++++++++++++++++               | 69% ~01m 46s      
   |+++++++++++++++++++++++++++++++++++               | 70% ~01m 42s      
   |++++++++++++++++++++++++++++++++++++              | 71% ~01m 39s      
   |++++++++++++++++++++++++++++++++++++              | 72% ~01m 35s      
   |+++++++++++++++++++++++++++++++++++++             | 73% ~01m 32s      
   |+++++++++++++++++++++++++++++++++++++             | 74% ~01m 29s      
   |++++++++++++++++++++++++++++++++++++++            | 75% ~01m 25s      
   |++++++++++++++++++++++++++++++++++++++            | 76% ~01m 22s      
   |+++++++++++++++++++++++++++++++++++++++           | 77% ~01m 18s      
   |+++++++++++++++++++++++++++++++++++++++          | 78% ~01m 15s      
   |++++++++++++++++++++++++++++++++++++++++          | 79% ~01m 12s      
   |++++++++++++++++++++++++++++++++++++++++         | 80% ~01m 08s      
   |+++++++++++++++++++++++++++++++++++++++++         | 81% ~01m 05s      
   |+++++++++++++++++++++++++++++++++++++++++         | 82% ~01m 01s      
   |++++++++++++++++++++++++++++++++++++++++++        | 83% ~58s          
   |++++++++++++++++++++++++++++++++++++++++++        | 84% ~54s          
   |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~51s          
   |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~48s          
   |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~44s          
   |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~41s          
   |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~37s          
   |+++++++++++++++++++++++++++++++++++++++++++++    | 90% ~34s          
   |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~31s          
   |++++++++++++++++++++++++++++++++++++++++++++++   | 92% ~27s          
   |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~24s          
   |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~20s          
   |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~17s          
   |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~14s          
   |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~10s          
   |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~07s          
   |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~03s          
   |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed = 05m 40s

DE_E16_ciliated_genotype

write.table(DE_E16_ciliated_genotype,"DE_E16_ciliated_genotype.txt",sep="\t")

DE_E16_ciliated_genotype$gene<-rownames(DE_E16_ciliated_genotype)
E16_ciliated_automatic_geneList<-DE_E16_ciliated_genotype$gene[DE_E16_ciliated_genotype$p_val_adj<0.001 & abs(DE_E16_ciliated_genotype$avg_logFC)>0.5 & abs(DE_E16_ciliated_genotype$pct.1-DE_E16_ciliated_genotype$pct.2)>0.15]

library(ggrepel)

#DE_P4_secretory_genotype$sig<-DE_P4_secretory_genotype$p_val_adj<0.001
DE_E16_ciliated_genotype$threshold<- ifelse(DE_E16_ciliated_genotype$avg_logFC>0 & DE_E16_ciliated_genotype$p_val_adj<0.001, "wt_enrich",ifelse(DE_E16_ciliated_genotype$avg_logFC<0 & DE_E16_ciliated_genotype$p_val_adj<0.001, "mut_enrich","NotSignificant" ) )
ggplot(DE_E16_ciliated_genotype, aes(avg_logFC, -log10(p_val_adj))) + #volcanoplot with avg_logFC versus p_val_adj
    geom_point(aes(col=threshold),size=0.2) + #add points colored by significance
  scale_color_manual(values=c("green", "black","magenta"))+
    ggtitle("E16Ciliated_wt/mut") + geom_text_repel(data=DE_E16_ciliated_genotype[DE_E16_ciliated_genotype$gene %in% E16_ciliated_automatic_geneList,], aes(label=gene), point.padding = 1, box.padding = .3) +
  labs(y = expression(-log[10]*" "*"adjusted pvalue"), x = "avg log fold change") + 
  theme(legend.title = element_blank(), legend.position = "top") 

PCD genes are not significantly differentially expressed between genotypes:

#DE_P4_secretory_genotype$sig<-DE_P4_secretory_genotype$p_val_adj<0.001
DE_E16_ciliated_genotype$threshold<- ifelse(DE_E16_ciliated_genotype$avg_logFC>0 & DE_E16_ciliated_genotype$p_val_adj<0.001, "wt_enrich",ifelse(DE_E16_ciliated_genotype$avg_logFC<0 & DE_E16_ciliated_genotype$p_val_adj<0.001, "mut_enrich","NotSignificant" ) )
ggplot(DE_E16_ciliated_genotype, aes(avg_logFC, -log10(p_val_adj))) + #volcanoplot with avg_logFC versus p_val_adj
    geom_point(aes(col=threshold),size=0.2) + #add points colored by significance
  scale_color_manual(values=c("green", "black","magenta"))+
    ggtitle("E16Ciliated_wt/mut") + geom_text_repel(data=DE_E16_ciliated_genotype[DE_E16_ciliated_genotype$gene %in% geneList$Primary.ciliary.dyskinesia,], aes(label=gene), point.padding = 1, box.padding = .3) +
  labs(y = expression(-log[10]*" "*"adjusted pvalue"), x = "avg log fold change") + 
  theme(legend.title = element_blank(), legend.position = "top") 

df_E16_epi<-FetchData(E16_Dec7v3_epi,c("Spdef","Creb3l1","Scgb3a2","Scgb1a1","Krt4","Krt13","Foxa3","Aqp3","Aqp4","Aqp5","Gp2","Sostdc1","Smoc2","Krt14","Krt15","Krt5","Rac2","Clic3","res.1.2","res.1.4","genotype","seq_group","specific_type","cell_type","Defb1","Lyz2","Ltf","Sftpa1","Sftpd","Sftpb","Slpi","Lcn2","Pigr","Muc5b","Muc5ac","Chil4","Muc1","Muc2","Muc4","Muc16","Muc20","Lbp","Cd14","Tlr4","Tlr2","Myd88","Ticam1","Itln1","Lgals3","Reg3g","Nod1","Nod2","Ddx58","Ifih1","Dhx58","Ccl5","Cxcl10","Cxcl2","Cxcl1","Pf4","Cxcl12","Cxcl14","Cxcl15","Cxcl16","Cxcl17","Ccl2","Ccl7","Ccl17","Ccl20","Ccl21a","Ccl25","Ccl27a","Ccl28","Cx3cl1","Il10","Tnf","S100a8","S100a9","Il6","Il18","Il1b","Il1rl1","Ccl11","Ccl24","Il33","Il25","Tslp","F2rl1","Retnla","Alox15","Alox5","Gata2","Tgfb2","Tgfb1","Ormdl3","Ptges","Ptgds","Ptgs2","Hpgds","Tbxas1","Areg"))

Error in FetchData(E16_Dec7v3_epi, c("Spdef", "Creb3l1", "Scgb3a2", "Scgb1a1",  : 
  Error: Chil4 not found

MicrobialSensing:

for (i in c("Lbp","Cd14","Tlr4","Tlr2","Myd88","Ticam1","Itln1","Reg3g","Lgals3","Nod1","Nod2","Ddx58","Ifih1","Dhx58"))
{
pdf(file = paste("Manuscript/MicrobialSensing_genotype/E16/",i,".pdf", sep = ""), width = 6, height = 5)
print(ggplot(df_E16_epi,aes_string(x="genotype",y=i))+facet_grid(.~cell_type)+geom_dotplot(binaxis="y",aes(fill=genotype),binwidth=0.05,stackdir="center",position=position_dodge(0.8), dotsize=0.2)+stat_compare_means(comparisons = list(c("wt", "mut")),method="wilcox.test",size=4,label="p.adj")+ stat_summary(aes(color=genotype),fun.data=mean_sdl, fun.args = list(mult=1), 
                 geom="pointrange",position=position_dodge(0.7))+ theme(axis.text.x = element_text(angle = 45,hjust=1),strip.text.x = element_text(size = 9, colour = "black", angle = 0)))
dev.off()
}

antimicrobial effectors:

for (i in c("Muc1","Muc4","Muc16","Muc20","Muc5b","Muc5ac","Muc2","Defb1","Lyz2","Ltf","Sftpa1","Sftpd","Sftpb","Slpi","Lcn2","Pigr","Chil4"))
{
pdf(file = paste("Manuscript/Effectors_genotype/E16/",i,".pdf", sep = ""), width = 6, height = 5)
print(ggplot(df_E16_epi,aes_string(x="genotype",y=i))+facet_grid(.~cell_type)+geom_dotplot(binaxis="y",aes(fill=genotype),binwidth=0.05,stackdir="center",position=position_dodge(0.8), dotsize=0.2)+stat_compare_means(comparisons = list(c("wt", "mut")),method="wilcox.test",size=4,label="p.adj")+ stat_summary(aes(color=genotype),fun.data=mean_sdl, fun.args = list(mult=1), 
                 geom="pointrange",position=position_dodge(0.7))+ theme(axis.text.x = element_text(angle = 45,hjust=1),strip.text.x = element_text(size = 9, colour = "black", angle = 0)))
dev.off()
}

Error in FUN(X[[i]], ...) : object 'Chil4' not found

chemokines:

for (i in c("Ccl5","Cxcl10","Cxcl2","Cxcl1","Pf4","Cxcl12","Cxcl14","Cxcl15","Cxcl16","Cxcl17","Ccl2","Ccl7","Ccl17","Ccl20","Ccl21a","Ccl25","Ccl27a","Ccl28","Cx3cl1"))
{
pdf(file = paste("Manuscript/chemokines_genotype/E16/",i,".pdf", sep = ""), width = 6, height = 5)
print(ggplot(df_E16_epi,aes_string(x="genotype",y=i))+facet_grid(.~cell_type)+geom_dotplot(binaxis="y",aes(fill=genotype),binwidth=0.05,stackdir="center",position=position_dodge(0.8), dotsize=0.2)+stat_compare_means(comparisons = list(c("wt", "mut")),method="wilcox.test",size=4,label="p.adj")+ stat_summary(aes(color=genotype),fun.data=mean_sdl, fun.args = list(mult=1), 
                 geom="pointrange",position=position_dodge(0.7))+ theme(axis.text.x = element_text(angle = 45,hjust=1),strip.text.x = element_text(size = 9, colour = "black", angle = 0)))
dev.off()
}

Th2:

for (i in c("Il10","Tnf","S100a8","S100a9","Il6","Il18","Il1b","Il1rl1","Ccl11","Ccl24","Il33","Il25","Tslp","F2rl1","Retnla","Alox15","Alox5","Gata2","Tgfb2","Tgfb1","Ormdl3","Ptges","Ptgds","Ptgs2","Hpgds","Tbxas1","Areg"))
{
pdf(file = paste("Manuscript/Th2_genotype/E16/",i,".pdf", sep = ""), width = 6, height = 5)
print(ggplot(df_E16_epi,aes_string(x="genotype",y=i))+facet_grid(.~cell_type)+geom_dotplot(binaxis="y",aes(fill=genotype),binwidth=0.05,stackdir="center",position=position_dodge(0.8), dotsize=0.2)+stat_compare_means(comparisons = list(c("wt", "mut")),method="wilcox.test",size=4,label="p.adj")+ stat_summary(aes(color=genotype),fun.data=mean_sdl, fun.args = list(mult=1), 
                 geom="pointrange",position=position_dodge(0.7))+ theme(axis.text.x = element_text(angle = 45,hjust=1),strip.text.x = element_text(size = 9, colour = "black", angle = 0)))
dev.off()
}

explore genes correlated with basal–>secretory changes:

c3 is the basal/secretory population

E16_Dec7v3_epi <- SetAllIdent(object = E16_Dec7v3_epi, id = "res.1.4")
E16_Dec7v3_epi_subc3<-SubsetData(object=E16_Dec7v3_epi,ident.use=c(3))
table(E16_Dec7v3_epi_subc3@meta.data$res.1.4)

E16_Dec7v3_epi_subc3@meta.data<-E16_Dec7v3_epi_subc3@meta.data[,-which(names(E16_Dec7v3_epi_subc3@meta.data) %in% c("res.0.8", "res.1.2", "res.1.4", "res.1.6"))] #remove old metadata
E16_Dec7v3_epi_subc3 <- ScaleData(object = E16_Dec7v3_epi_subc3)

E16_Dec7v3_epi_subc3 <- FindVariableGenes(object = E16_Dec7v3_epi_subc3, do.plot = TRUE, x.low.cutoff=0.1,x.high.cutoff = Inf, y.cutoff = 0.5)

E16_Dec7v3_epi_subc3 <- RunPCA(object = E16_Dec7v3_epi_subc3, do.print = FALSE)
PCAPlot(E16_Dec7v3_epi_subc3)

PCElbowPlot(object = E16_Dec7v3_epi_subc3)

n.pcs.sub3 = 13

resolution parameter sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater number of clusters.

res.used <- 1.2

E16_Dec7v3_epi_subc3 <- FindClusters(object = E16_Dec7v3_epi_subc3, reduction.type = "pca", dims.use = 1:n.pcs.sub3, 
                     resolution = res.used, print.output = 0, force.recalc = T)

E16_Dec7v3_epi_subc3 <- RunTSNE(object = E16_Dec7v3_epi_subc3, dims.use = 1:n.pcs.sub3, perplexity=30)

TSNEPlot(object = E16_Dec7v3_epi_subc3, do.label = T)

res.used <- 0.8

TSNEPlot(object = E16_Dec7v3_epi_subc3, do.label = T)

res.used <- 0.6

TSNEPlot(object = E16_Dec7v3_epi_subc3, do.label = T)

E16_Dec7v3_epi_subc3 <- SetAllIdent(object = E16_Dec7v3_epi_subc3, id = "res.0.6")
DoHeatmap(object = E16_Dec7v3_epi_subc3, genes.use = c("Epcam","Krt8","Trp63","Krt5","Mki67","Top2a","Smoc2","Ccl20","Sostdc1","Bmp7","Clic3","Cldn10","Tspan33","Ehf","Sfta2","Crip2","Msln","Cyp2s1","Cldn3","Cldn7","Cldn4","AU021092","Tspan1","Chad","Tspan13","Klk10","Klk11","Klk13","Ces1d","Krt4","Krt13","Creb3l1"), 
    slim.col.label = TRUE, group.label.rot = TRUE,use.scaled = T
  )

cor(E16_Dec7v3_epi_subc3@scale.data["Trp63",], E16_Dec7v3_epi_subc3@scale.data["Creb3l1",])

cor(E16_Dec7v3_epi_subc3@data["Trp63",], E16_Dec7v3_epi_subc3@data["Spdef",])

Trp63_cor<-apply(E16_Dec7v3_epi_subc3@scale.data, 1, function(x) cor(E16_Dec7v3_epi_subc3@scale.data["Trp63",],x)) 

head(Trp63_cor) 

min(Trp63_cor,na.rm = T) 

Trp63_cor_order<-order(Trp63_cor,decreasing=T) 

head(Trp63_cor[Trp63_cor_order],20)

sum(is.na(Trp63_cor[Trp63_cor_order]))

Trp63_cor[Trp63_cor_order][15706:15686]

tail(Trp63_cor[Trp63_cor_order],20)

cor.test(E16_Dec7v3_epi_subc3@scale.data["Spdef",], E16_Dec7v3_epi_subc3@scale.data["Trp63",],method="pearson")

cor.test(E16_Dec7v3_epi_subc3@scale.data["Spdef",], E16_Dec7v3_epi_subc3@scale.data["Trp63",],method="kendall")

Trp63_cor_test<-apply(E16_Dec7v3_epi_subc3@scale.data, 1, function(x) cor.test(E16_Dec7v3_epi_subc3@scale.data["Trp63",],x,method="pearson")) 

df_Trp63_corTest<-as.data.frame(do.call(rbind, Trp63_cor_test))

df_order_Trp63_corTest<-df_Trp63_corTest[order(unlist(df_Trp63_corTest$estimate)),]

tidy_Trp63_cor<-cbind(df_order_Trp63_corTest$estimate,df_order_Trp63_corTest$p.value)

colnames(tidy_Trp63_cor)<-c("cor","p.value")

head(tidy_Trp63_cor,20)

tidy_Trp63_cor[15706:15686,]

tail(tidy_Trp63_cor)

Krt8_cor_test<-apply(E16_Dec7v3_epi_subc3@scale.data, 1, function(x) cor.test(E16_Dec7v3_epi_subc3@scale.data["Krt8",],x,method="pearson")) 

df_Krt8_corTest<-as.data.frame(do.call(rbind, Krt8_cor_test))

df_order_Krt8_corTest<-df_Krt8_corTest[order(unlist(df_Krt8_corTest$estimate)),]

sum(is.na(df_order_Krt8_corTest$estimate))

tidy_Krt8_cor<-cbind(df_order_Krt8_corTest$estimate,df_order_Krt8_corTest$p.value)

colnames(tidy_Krt8_cor)<-c("cor","p.value")

head(tidy_Krt8_cor,20)

tidy_Krt8_cor[15666:15706,]

cor.test(E16_Dec7v3_epi_subc3@scale.data["Krt8",E16_Dec7v3_epi_subc3@meta.data$genotype=="wt"], E16_Dec7v3_epi_subc3@scale.data["Numb",E16_Dec7v3_epi_subc3@meta.data$genotype=="wt"],method="pearson")

cor.test(E16_Dec7v3_epi_subc3@scale.data["Krt8",E16_Dec7v3_epi_subc3@meta.data$genotype=="mut"], E16_Dec7v3_epi_subc3@scale.data["Numb",E16_Dec7v3_epi_subc3@meta.data$genotype=="mut"],method="pearson")